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ABSTRACT 


A  boundary  integral  solution  to  the  two-dimensional,  free-surface  water  wave  problem  is 
presented.  The  mathematical  formulation  involves  application  of  potential  theory  and 
appropriate  initial  and  boundary  conditions  to  resolve  the  progression  of  the  linear  free-surface 
waves.  The  solution  to  the  potential  flow  problem  is  represented,  through  Green's  theorem,  by  a 
boundary  integral  method  that  is  approximated  via  linear  boundary  panels.  The  resulting  system 
of  algebraic  equations  is  solved  for  required  flow  parameters.  The  free  surface  is  tracked  at  each 
time  level  by  numerical  integration  of  the  linearized  free-surface  boundary  conditions.  Despite 
the  use  of  linear  panels,  numerical  results  compare  favorably  with  the  exact  linear  theory  and 
indicate  that  the  computational  scheme  is  able  to  track  the  development  and  propagation  of  steep 
waves  over  long  periods  of  time.  An  analysis  is  also  provided  to  aid  in  the  development  and 
numerical  implementation  of  artificial  or  perfectly  absorbing  boundary  conditions  at  the  vertical, 
imaginary,  truncated  boundaries. 
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A  BOUNDARY  INTEGRAL  METHOD 
FOR  TWO-DIMENSIONAL  WATER  WAVES 


INTRODUCTION 

Recent  experiments  by  Kuria  et  al.,1  as  well  as  earlier  research  by  other  investigators,2'9 
have  shown  that  water  wave  interaction  can  have  a  significant  hydrodynamic  effect  on  marine 
structures.  Bodies  of  water  always  have  surface  waves,  which  are  due  to  the  effects  of  forces 
acting  on  the  fluid  and  trying  to  deform  it.  The  size  and  form  of  a  surface  wave  will  depend  on  the 
size  and  magnitude  of  the  wave-generating  mechanism,  while  its  propagation  will  be  influenced  by 
gravity  and  surface  tension. 

The  theory  of  water  waves  is  provided  in  great  detail  by  Dean  and  Dalrymple.2  The 
analytical  solution  to  the  linear  potential  flow  equations  assumes  that  the  vertical  lateral  boundary 
conditions  are  periodic  in  space  and  time.  While  it  is  fairly  easy  to  implement  the  spatial  periodic 
boundary  condition,  its  quite  difficult  to  numerically  implement  the  temporal  periodicity  condition. 
Longuet-Higgins  and  Cokelet3  used  coordinate  transformation  to  convert  the  computational  domain 
to  a  simple  closed  contour.  Vinje  and  Brevig4  have  solved  the  nonlinear  two-dimensional  wave 
problem  in  physical  space.  To  satisfy  periodic  conditions  at  the  truncated,  lateral  boundaries,  they 
imposed  periodicity  in  both  velocity  potential  and  stream  function.  A  more  general  boundary 
integral  solution  to  the  nonlinear  two-dimensional  wave  equation  was  provided  by  Sen  et  al.5  In 
this  work,  the  authors  used  linear  panels  to  simulate  motions  of  two-dimensional  large  floating 
bodies  in  waves. 

To  avoid  the  problems  associated  with  lack  of  continuity  at  linear  panel  boundaries,  Sen6 
developed  a  cubic-spline  boundary  integral  method.  With  a  cubic  representation  of  the  velocity 
potential  in  a  panel,  it  is  possible  to  enforce  continuous  velocity  at  the  panel  boundaries  and  comer 
points.  Xii7  used  the  boundary  integral  method  to  obtain  a  numerical  solution  to  the  fully  nonlinear 
three-dimensional  water  wave  problem  in  a  tank.  This  three-dimensional  boundary  element 
method  is  based  on  biquadratic  isoparametric  curvilinear  elements  with  an  efficient  adaptive 
numerical  quadrature  scheme.  After  every  three  or  four  time  steps,  Chebyshev  smoothing  was 
used  in  alternating  directions  to  remove  the  three-dimensional  wave  instabilities. 

In  the  research  documented  here,  the  boundary  integral  method  is  used  to  simulate  the 
linear,  two-dimensional  wave  equation  in  the  time  domain.  This  computational  formulation  is 
based  on  Green’s  formula  for  harmonic  functions  applied  to  a  truncated  fluid  domain.  To  obtain  a 
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system  of  algebraic  equations  from  the  Fredholm  integral  equations,  the  proposed  approach  applies 
the  zeroth-order  linear  panels  along  the  control  boundaries.  The  resulting  system  of  equations  is 
then  solved  for  the  normal  velocity  and  velocity  potential.  The  tangential  velocity  is  obtained  by 
numerical  differentiation  of  the  velocity  potential.  The  free  surface  is  advanced  in  time  by 
numerical  integration  of  the  linearized  free-surface  kinematic  and  dynamic  boundary  conditions. 

To  start  the  solution,  an  excitation  Airy  wave  velocity  potential  is  imposed  on  the  leading  upstream 
vertical  boundary.  At  the  downstream  truncation  boundary,  either  the  spatial  and  temporal 
periodicity  condition  or  the  Sommerfeld  radiation  boundary  condition  is  applied. 
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MATHEMATICAL  FORMULATION 


The  theory  and  governing  physics  of  gravity  water  waves  and  ffee-surface  flow  are 
presented  in  Dean  and  Dalrymple.2  In  the  physical  geometry  defined  in  figure  1,  boundary  dDCi  is 
the  location  of  the  wave  generator  and  dDc  is  a  fictitious  boundary  placed  there  to  make  the 
computational  domain  finite.  Boundary  dDF  is  the  free  surface  and  dDD  represents  a  horizontal, 
nonporous  solid  boundary. 

The  physical  coordinates  (x,y)  represent  a  point  with  a  corresponding  velocity,  (w,v). 

The  velocity  components  are  related  to  the  velocity  potential  <p  through 


x=0  dDD  x=L 


Figure  1.  Geometric  Definitions 
We  give  the  following  definitions  of  the  domain  D  and  its  boundary  dD : 
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(2) 


D  =  \{x,y)\0<x<L,-d<y<T](x,t )}  , 
dDF  =  {(x,y)|0<x<i;,y  =  T7(x,r)}, 

dDCi  =  {(*,y)  | x  =  L,-d<y<  Ti(x,t)}  , 

9Dd  =  {0,y)  |0<x  <L,y  =  ~d\  , 

dDCi  =  {(x,y)|x  =  0,-</<  y<ri(x,t)}. 


With  the  assumption  of  irrotational  motion  and  an  incompressible  fluid,  a  velocity  potential  exists 
that  must  satisfy  the  governing  equation  of  mass  conservation.  In  terms  of  the  velocity  potential, 
the  continuity  equation  leads  to  the  Laplace  equation,  which  must  hold  throughout  the  fluid: 

V->  =  0  in  D.  (3) 

The  free-surface  kinematic  boundary  condition  requires  that  a  fluid  particle  on  the  surface 
remain  on  the  free  surface.  That  is, 

~{y-T](x,t)}  =  0  on  y(t)  =  rj(x,t)  ,  (4) 


which,  on  expanding,  gives 


chi  _d(j)  d(f)  dri 
dt  dy  dx  dx 


y(t)  =  rj(x,t). 


(5) 


The  pressure  within  the  fluid  domain  D  and  on  all  the  boundaries  must  conform  to  Bernoulli’s 
equation.  The  Bernoulli  equation  is  therefore  used  to  obtain  the  dynamic  free-surface  boundary 
condition.  In  its  irrotational  form,  this  boundary  condition  is 


d(f> 

dt 


1 

2 


+ 


A, 


~gl  on  y(t)  =  Tj(x,t). 


(6) 


On  the  rigid,  nonporous  bottom  boundary,  there  can  be  no  fluid  velocity  component  normal 
to  the  boundary.  The  flow  must  be  entirely  tangential  to  the  boundary,  so  that  the  bottom  surface 
is  always  a  streamline.  This  can  be  mathematically  expressed  as 
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^  =  0  on  3D, 
an 


(7) 


D  > 


where  n  is  the  unit  normal  pointing  outward  from  the  computational  domain. 
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NUMERICAL  ANALYSIS 


Solutions  to  the  harmonic  velocity  potential  may  be  represented  by  its  boundary  data. 

In  a  two-dimensional  region,  Green’s  second  identity  gives 

«x)=7rf-f  «©|-G(x.4)-G(x,y|;<t©  ds&  ,  (8) 

Q(x)JdL  dn  dn  J 


where  d / dn  denotes  the  derivative  along  the  outward  facing  normal  n  on  9D,  and 
G(x,  $,)  =  In  is  the  kernel  Green's  function,  where 

R(x,0  =  |x  -  5|  =  [(*  -  if  +  (y  -  r]f]K  (9) 

is  the  distance  between  the  points  x  and  c,  on  the  boundary.  The  expression x  =  xi  +  yj  defines 
any  point  in  D  =  D  +  dD ,  \  +  rij  defines  any  point  on  dD,  and 

k,  x  g  smooth  part  of  dD 

Q(x)  =  ■  internal  angle,  xe  comer  of  dD  .  (10) 

2k,  x  e  D 

v 


In  boundary  integral  calculations,  the  computational  domain  is  decomposed  into  many 
panels  and  an  interpolation  polynomial  is  used  to  approximate  the  velocity  potential  and  normal 
velocity  on  each  panel.  The  accuracy  of  the  computational  technique  is  improved  either  by  an 
increase  in  the  number  of  computational  panels  or  by  the  order  of  the  interpolation  polynomial.  If 
N  is  the  order  of  the  interpolation  polynomial  and  iVy(q)  represent  the  cardinal  functions,  then  the 
velocity  potential  at  any  point  on  the  computational  domain  is  approximated  by 


and 


0©=5X©9>; 


vft  pi 


(11) 


where  s  indicates  the  particular  time  level,  i.e.,  t  -  s(Al)  (At  is  the  time  step  and  Af(^)  are 
appropriate  basis  functions).  When  linear  panels  are  used,  0  and  d(f)/ dn  are  ultimately  assumed 
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to  be  constant  in  each  panel.  In  this  work,  it  will  therefore  be  assumed  that  Af(£)  are  always 
constant  over  every  boundary  panel.  When  equation  (1 1)  is  substituted  into  equation  (8),  the  finite 
dimensional  approximation  to  0(x)  is  found  to  be 


= 7J-7 J  (z nj® <e) 1 4 1 ^ 1 & - c(*. ■ d S *0© 

i2(X),9£>  v>  1  y10”  V/=l 


(12) 


j=l  VoD 


—if  f  <*& 


Q(x) 


If  we  restrict  xedD  and  define  a  projection  (collocation)  operator,  PN,  such  that 

PN(t>(x)  =  <t>(xk)  , 

the  following  approximation  is  obtained: 


^«v«=  f,^)  <?■ = = < 

i=i  y=i 


where  x*  e  oD  are  the  boundary  collocation  points  and 


(13) 


)=<?,*,  1  <(j,k)<N. 


(14) 


Applying  the  definition  of  PN  to  equation  (12),  we  obtain 


PM*)  =  XPJ  7^  1  NJ&-?-G(x&)ds(g) 

j=  i  y^Ax) 


<P 


'  3D 


rj  ■ 


(15) 


To  be  able  to  evaluate  the  integrals  in  equation  (15),  we  partition  dD  into  N  panels  (line 
segments)  such  that 


dD  =  \JdD, 

/= 1 


(16) 
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and 


AT,®  = 


1,  i,sSDj 

0,  otherwise. 


(17) 


The  integrals  in  equation  (15)  can  therefore  be  written  as 


f Nj(%)^-G{xfyds^)  =  X  f N/?;)~G(x,t;)ds(Z) 
L  d*  tT  &'  dn 


=  f  |-G(x^)^)5 


(18) 


and,  likewise, 


|G(x,|)A',®*«)=  Jg(x,®<6®  .  (19) 

5D  dDj 


The  expressions  on  the  right-hand  side  of  equation  (15)  may  now  be  evaluated  as 


J-  f  N 0^-G(x&)ds<S) )  =  -4-  J  2-G(x,y  <6® 

£Xx)i  *  I  nrxji  * 


y  Q(x\)/Didn 
=  1Hki 

K  *J 


X=Xk 


(20) 


and 


where 


^  J Gd&Jf/S) *©  =  -4-  J Gd,5)*®| 


Q(x) 


(21) 


Gt;=  jG(x,§)<6(«|„.  , 

dDj 

*„ = j  u 

dD  , 


(22) 
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and  Q(x*J  =  n .  For  our  purposes,  the  collocation  (control)  points  x'k  will  be  chosen  to  be 
midpoints  of  dDk . 

Upon  substituting  equations  (13),  (20)  and  (21)  into  equation  (15),  the  following  system  of 
equations  is  obtained: 


'Pj'-iXv',’.  l£k<N  ,  (23) 

J=  i  M 

or,  equivalently. 


siJv;='Lniyi-'Loi/¥-,  i <k<N 


(24) 


j= 1 


J=  1 


which  may  be  written  as 

f,Hv<P'-f1Gi/rj=0,  1<  k<N,  (25) 

J=  1  >1 

where 


Hkj  —  Hkj  7c8kj  . 


(26) 


We  next  particularize  the  system  of  equations  relative  to  two-dimensional  wave  generation  in  a 
rectangular  cavity.  The  domain  is  given  in  figure  1.  We  partition  the  respective  parts  of  dD  as 
follows: 


bdf 

c1 

VI 

VI 

<p  -  specified , 

np\  + 1  <  j  <  np2 

0  -  specified , 

dDD 

np2  + 1  <  j  <  np3 

^  -  specified , 
an 

dDCl 

np3  + 1  <  j  <  N 

0  -  specified , 
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where  np\,  np2  -  np\,  np3  -  np2,  and  N  -  np3  are,  respectively,  the  number  of  boundary  panels  on 
dD  f ,  dDq  ,  dD  D ,  and  dDc^ . 


Substituting  the  governing  information  given  in  equation  (27)  into  equation  (25),  including 
the  known  boundary  conditions,  results  in  the  following  system  of  algebraic  equations: 


np3  np2  N 

%GkJxif) 

J=np2+ 1  j=- 1  j^np3+ 1 


(28) 


npl 


np3 


= is*^ 

j=l  j=npl+l  J=np2+\ 


which  may  be  written  symbolically  as 


Axs  =  b  , 


(29) 


where  the  elements  of  xs  are  given  by 


Vi 


x 


s 


(30) 


Vn 


The  elements  of  A  are 
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A  = 


H \jnp2+\ 

- 

H  np2,np2+\ 

^  np2pp3 

£> 

11  np2+\/ip2+\ 

£> 

***  ^  np2+\,np3 

np3,np2+\ 

^  np3,np3 

ZJ 

r\p3+\,np2+\ 

^  np3+\pp3 

^N^p2+\ 

^N/ip3 

^l,np3+l 


'^Npp2+\  ""  &N,N 


(31) 


and  the  elements  of  b  are 


np2  np3  N 

h=-'L«vvJ+  ?,Gvr,-  ak<N 

j=  1  j=np2+\  j=np2+l 


(32) 


We  note  that  the  sub-block 


H  H 

**  np2+  \,np2+ 1  *  *  *  11  np2+\,np3 


H 


np3,np2+l 


H 


np3tnp3 


(33) 


is  the  only  portion  of  A  that  has  K  subtracted  along  its  main  diagonal. 

From  Nachbin,8  the  matrix,  as  defined  by  equation  (26),  must  have  the  following  property* 
for  any  fixed  fc  — »  i: 


k,  -  = -S  • 


M 

j*i 


j*i 


(34) 


Therefore, 


*The  code  we  used  in  our  simulations  was  modified  to  reflect  this  property. 
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(35) 


N 

H,i=-'LHtj+x  ■ 

M 

j*i 

If  equation  (25)  is  to  model  a  constant  potential  (f>  in  each  panel,  then  we  must  have  y/j  =0,  and 
therefore 


p;=0  (36) 

M 

for  any  constant  vector  .  Equations  (34)  and  (35)  follow  immediately  from  equation  (36). 
Equation  (35)  provides  a  convenient  way  of  evaluating  the  second  part  of  equation  (22)  when 
k  =  j. 


In  order  to  solve  the  system  of  algebraic  equations  in  equation  (29),  expressions  are 
required  for  the  integrals  in  equation  (22).  To  maintain  legibility,  the  following  function 
definitions  are  introduced: 


and 


¥(x)  =  |  In  Rdsfc) 

*>J 


(37) 


<fc(x)=  |  |-(lni?)<fc(S)  • 


(38) 


From  figure  2,  it  can  be  seen  that  on  dDj ,  r  =  x  -  £  =  (x  -  £)i  +  (y  -  rj)  j ,  and  therefore 
R=(r-r/>  =[(x-£)2+(y-T1)2/\ 

Using  figure  2,  we  may  next  rewrite  the  integrand  in  equation  (38)  as 

4<M)  =  n-V{(ln*)=-5-V{i  =  — .  (39) 

an  K  R 

Equation  (38)  thus  becomes 
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(40) 


®00  =  -  /  • 

4  R 

Equation  (22)  may  then  be  expressed  in  terms  of  equations  (37)  and  (40)  as 


GkJ  =  'F(xl) 

and  (41) 


The  integrals  in  equations  (37)  and  (40)  are  evaluated  using  the  ideas  introduced  by  Faltinsen.9 
The  development  that  follows  is  based  on  the  geometry  of  figure  2.  The  subscripts  L  and  R  , 
respectively,  denote  the  left  and  right  endpoints  of  the  boundary  element  dD}\  s  =  s(£,)  and 
n  =  n(q) ,  respectively,  denote  the  unit  tangent  and  normal  to  dD} ,  where  s  points  from  c,L  to  t,R 
and  n  =  k  a  s  points  outward  to  dDj  at  (£77) ;  and  (—s,—  n)  are  the  components  of  r  relative  to 
the  basis  (s,n) .  In  defining  the  outward  normal  n,  it  is  assumed  that  the  boundary  element  dDJ  is 
oriented  such  that,  relative  to  the  direction  s,  the  interior  of  the  region  D  is  on  the  right  (see  figure 
2).  The  following  relationships  are  then  obtained: 


st<s<sR  , 


.  ,  A77  Art 

s  =  — n - j,  11  = - n - J  , 

As  As  As  As 


(AZ,At})  =  ((£r-Zl),(tir-til)),  As=^/(A£)2  +  (A77)2 


n  =  -n  r  =  -n  rJ=-n-r1,  s  =  -s  r 


r  =  -«n-ss,  and  R  =  n2  +  s2 


(42a) 

(42b) 

(42c) 

(42d) 

(42e) 


Equation  (37)  therefore  becomes 
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where  sL=- s  •  rL  and  sR  =  -s  -  rR .  The  normals  may  be  defined  either  by  n  = 
n  =  -  n  rR.  Similarly,  equation  (40)  becomes 


*R 


*r 


®W=-I^r^5)  =  I 


n 2  +  52 


i/5, 


tan  1 

(s  ^ 
iJL 

-tan  1 

(s  ) 
1L 

x  €  dDj  , 


(44) 


where  the  quantities  sL,  sR,  and  n  all  depend  on  x  through  rL  and  rR.  Of  particular  interest  is  the 
evaluation  of  equation  (40)  when  x  =  x* ,  x*  e  dD  . ,  and  x*  ^  or  .  One  obtains  different 
values  for  O(x)^. ,  depending  on  how  the  limit  x  -»  x”  is  taken.  If  the  point  x*  is  approached 
through  positive  values  of  n,  then  for  values  of  X  sufficiently  close  to  x*,  it  is  easy  to  see 
(since sL  <  0  and  sR  >  0)  that 


d>(x)  =  5T.  (45) 

X— »X* 

h0+ 


If  the  point  x*  is  approached  through  negative  values  of  n ,  then  for  values  of  x  sufficiently  close 
to  x*, 


O(x)  =  -it .  (46) 

x->x* 
n->  0“ 


The  limit  taken  through  positive  values  of  n  means  that  x  — >x”  from  points  within  D,  and, 
likewise,  the  limit  taken  through  negative  values  of  n  means  that  x  -»  x*  from  points  outside  D. 
Finally,  if  the  point  x*  is  approached  through  points  on  dD} ,  then  n  =  0 ,  and  since  sL  and  sR  are 
independent  of  n,  we  have 


O(x)  =0.  (47) 

x-*x'jjtedDj 

ihO 


It  is  therefore  appropriate  to  assume  that  when  equation  (8)  is  used  for  the  potential  <j>,  the 
latter  of  these  three  approaches  is  the  correct  way  to  evaluate  equation  (40)  when  x  =  x* , 
x*  e  dDj ,  and  x*  or  l)R .  This  method  of  evaluating  the  limits  is  consistent  with  the  finite 
dimensional  approximation  to  the  flow  parameters.  In  these  calculations,  X  has  been  explicitly 
restricted  to  points  on  dD .  The  proposed  method  of  evaluating  the  limit  is  consistent  with 
equation  (36).  When  the  control  points  are  chosen  to  be  the  midpoints  of  dDp  it  can  be  shown  that 
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for  each  fixed  value  of  i,  1  <  i  <  N .  Similar  results  were  also  obtained  by  Nachbin8  while  using  a 
conjugate  harmonic  representation  for  the  integrand  presented  in  equation  (40).  Results  from 
boundary  integral  methods  are  known  to  be  independent  of  the  choice  of  boundary  collocation 
points  or  boundary  element  basis. 
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TEST  PROBLEM  ALGORITHM 


The  accuracy  of  the  proposed  method  and  Green  function  integrals  is  tested  with  a  simple 
problem,  which  is  that  of  potential  flow  in  a  fixed,  rectangular  computational  domain.  The 
governing  equation  is  the  Laplace  equation  on  a  (0,l)x  (-1,0)  rectangular  geometry  with  the 
following  boundary  conditions: 


0  =  2  on  dDF  -  {(*,>/)  |  y  =  o}  , 
0  =  2  on  dDCi  =  {(x,j)  |  x  =  lj  , 

~  =  0  on  9Dd=  {(*,>>)  |  y  =  -l}  , 
0  =  2  on  dDCi  =  |  x  =  oj 


(49) 


The  unique  solution  to  this  problem  is  the  constant  potential  0(x,  y,t)  =  2.  To  test  the  applicability 
of  the  above  integral  equations,  the  problem  is  solved  with  10  equally  spaced  panels  on  each  of  the 
four  boundaries.  Results  are  presented  in  figure  3  for  control  points  located  at  a  quarter,  mid,  and 
three  quarters  of  the  panel.  As  a  measure  of  the  accuracy  of  the  numerical  solution,  we  have 
graphed  d(j>/dy\  f  at  the  collocation  (control)  points  for  different  partitioning  of  dD.  The  exact 
value  is  d(j)/ dy\  _o  =  0,  Vx,  which  is  in  excellent  agreement  with  the  values  shown  in  figure  3, 
indicating  that  the  results  obtained  from  these  integrals  are  independent  of  the  location  of  the 
control  points. 

Another  issue  to  be  discussed  is  the  accuracy  of  the  numerical  method  when  two  adjacent 
panels  are  of  different  lengths.  It  should  be  noted  that  the  accuracy  of  the  linear  panel  numerical 
method  is  usually  improved  by  increasing  the  number  of  panels.  The  question  is  how  to  arrange 
the  same  number  of  panels  in  the  computational  domain  so  that  their  physical  geometry  is 
adequately  defined.  If  the  physical  domain  is  properly  defined,  the  length  of  adjacent  panels  will 
have  negligible  effect  on  the  results.  In  figure  4,  the  Chebyshev-Gauss-Lobatto  points  were  used 
to  define  the  physical  geometry.  These  points  are  known  to  produce  a  grid  that  is  coarse  at  the  end 
points  and  sparse  at  the  middle.  The  control  points  were  placed  at  the  middle  of  each  panel  in  each 
case  presented.  Figure  4  shows  that  despite  the  fact  that  the  panels  were  of  different  lengths,  the 
results  completely  agree  with  those  for  the  equally  spaced  panels  of  figure  3.  This  indicates  that 
the  length  of  adjacent  panels  does  not  significantly  affect  the  accuracy  of  the  computational  domain. 
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Figure  4.  Computed  Values  of  dfy/dy  with  Different  Panel  Lengths 
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NUMERICAL  RESULTS 


Next  we  consider  the  numerical  solution  to  the  linearized  water  wave  problem  as  defined  by 
the  Fredholm  integral  relation  presented  in  equation  (8).  To  simulate  the  propagation  of  an  Airy 
wave  in  the  computational  domain,  the  initial  value  of  the  velocity  potential  on  the  undisturbed  free 
surface,  obtained  from  the  periodic  linear  solution,  is 


m= 


HX  cosh(&0+  d)) 
IT  sinh(fcf) 


sin(Lc-crr)  . 


(50) 


This  velocity  potential  excitation  is  chosen  to  correspond  to  an  Airy  wave  of  height  H,  wavelength 
X,  and  period  T  propagating  in  the  positive  x-direction.  Other  parameters  are  k  =  2 tt/A,  which 
represents  the  wavenumber,  and  <7  =  2 k/T,  which  is  the  wave  frequency. 

To  complete  specifications  of  the  boundary  value  problem,  conditions  must  also  be 
specified  on  the  remaining  lateral  boundaries,  dDc^  and  dDc  .  For  the  initial  linear  solution,  the 
Airy  wave  excitation  described  by  equation  (50)  is  applied  on  the  upstream  vertical  boundary,  and 
spatial  and  temporal  periodicity  is  maintained  on  the  downstream  boundary.  For  the  wave 
excitation  described  by  equation  (50),  periodicity  is  achieved  by  applying  the  same  potential 
downstream. 

Numerically,  the  free-surface  profile  and  velocity  potential  must  be  recovered  after  every 
period.  The  governing  equation  and  boundary  conditions  are  to  be  evaluated  using  the  boundary 
integral  methods  presented  in  this  memorandum.  To  advance  the  solution  in  time,  the  fourth-order 
Adams-Bashforth  (AB4)  method  was  used.  The  fourth-order  Runge-Kutta  (RK4)  method  was 
also  tested  and  compared  to  AB4;  both  methods  were  shown  to  give  the  same  results.  Since  AB4 
only  needs  calculations  from  previous  time  levels  and,  unlike  RK4,  does  not  require  solving 
Laplace’s  equation  four  times  at  every  time- level,  AB4  was  found  to  be  more  computationally 
efficient  and  was  thus  used  for  the  numerical  calculations  presented  in  this  section.  To  start  the 
solution,  the  lower  order  Adams-Bashforth  method  was  used  with  correction  until  there  was 
enough  information  to  use  AB4  without  correcting. 

In  water  wave  theory,  the  free-surface  boundary  conditions  must  be  satisfied  on  the  free 
surface  y  =  ri(x,t),  which  is  a  priori  unknown.  A  convenient  method  for  evaluating  these 
conditions  is  to  expand  them  at  y  =  0  (a  known  location)  by  the  truncated  Taylor  series.  The 
dynamic  free-surface  boundary  condition,  equation  (6),  therefore  becomes 
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(51) 


/ 

V 


^  +  1 
,  dt  2 


+ 


=  0, 


where  the  pressure  is  taken  to  be  ambient  on  the  free  surface.  Retaining  the  terms  that  are  linear  in 
the  small  parameter  T|,  and  recalling  that  r\  is  only  a  function  of  x  and  t,  the  linearized  dynamic  free- 
surface  boundary  condition  becomes 


f 

V 


d(f) 


) 


y=0 


=  0  . 


The  kinetic  free-surface  boundary  condition  is  similarly  linearized  to  become 


drj  __  d(j) 

*  b  y=0 


(52) 


(53) 


Since  the  linearized  equations  are  applied  on  the  undisturbed  free  surface  y  =  0,  the  entire 
computational  boundary  is  independent  of  time.  Consequently,  the  influence  coefficients  in  the 
resulting  system  of  linear  equations  remain  unchanged  in  time. 

Computed  free-surface  progressions  with  time  for  a  wave  with  an  amplitude  of  0.4  m, 
wavelength  of  2.46  m  (steepness  of  0.33),  depth  of  4.0  m,  and  overall  axial  length  of  4 
wavelengths  are  presented  in  figure  5.  A  time  step,  At  of  0.001  s,  was  used  in  these  simulations. 
The  other  parameters,  such  as  wave  frequency  and  period,  were  obtained  from  the  linear  water 
wave  theory.  The  results  show  that  the  proposed  method  can  obtain  wave  propagation  over  a  long 
period  of  time. 

Figure  6  shows  results  for  free-surface  progressions  near  the  upstream  boundary  at  x  = 
0.44  m,  near  the  center  at  x  =  4.87  m,  and  near  the  downstream  boundary  at  x  =  9.4  m.  Even  near 
the  boundaries,  it  is  seen  that  the  computed  results  compare  favorably  with  the  exact  linear 
solution. 
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The  challenge  of  the  water  wave  problem  is  to  be  able  to  simulate  the  fully  nonlinear  wave. 
The  linear  study  presented  in  this  memorandum  involves  the  first  phase  of  development,  which  is 
focused  on  determining  the  appropriate  boundary  conditions.  Though  the  free-surface  boundary 
conditions  are  well  known,  the  free  surface  itself  is  not  defined.  Another  boundary  problem 
concerns  the  truncated  upstream  and  downstream  boundaries.  Though  these  boundaries  are 
known,  their  boundary  conditions  still  remain  a  major  area  of  research.  The  actual  problem  is,  in 
fact,  due  to  these  boundaries  being  artificial  rather  than  real;  they  are  there  only  to  truncate  the 
computational  domain. 

In  this  work,  the  upstream  boundary  {at  x  =  0)  is  simulated  by  a  piston  wavemaker.  A 
velocity  potential  is  therefore  applied  to  this  boundary  at  all  times.  Although  another  alternative 
would  be  to  impose  the  wavemaker  velocity  on  this  boundary,  the  velocity  potential  alternative  is 
preferred  here  because  application  of  the  normal  velocity  would  require  a  moving  upstream 
boundary.  A  stationary  upstream  boundary  allows  one  to  use  a  fully  Eulerian  computational 
scheme. 

For  the  results  presented  in  figures  5  and  6,  the  linear  theory  was  used  to  obtain  the 
downstream  boundary  condition.  While  this  works  well  for  linear  waves,  convergence  problems 
would  be  expected  with  nonlinear  waves.  To  achieve  a  more  open  downstream  boundary 
condition,  the  Sommerfeld  radiation  condition  is  proposed.  The  radiation  condition  is  derived 
from  the  wave  equation 


dt2  dx2  ’ 


(54) 


where  0  is  any  quantity  and  C  is  the  phase  velocity  of  the  waves.  The  above  equation  can  be 
written  as 


( 

V 


£ 
dx  J 


0  =  0  . 


(55) 


The  first  term  on  the  left-hand  side  of  equation  (55)  represents  a  wave  traveling  in  the  positive  x- 
direction  with  a  speed  of  C.  For  this  wave,  the  above  equation  is  solved  to  give  the  radiation 
downstream  boundary  condition  as  follows: 
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(56) 


d(f)  _  ^d0 
dt  dx 

Based  on  the  same  parameters  as  those  of  figure  5,  the  downstream  results  from  this  boundary 
condition,  as  presented  in  figure  7,  compare  favorably  to  the  free-surface  progressions  of  figure  5. 

The  radiation  boundary  condition  has  been  tested  with  waves  of  different  depth  and 
steepness.  Results  for  the  wave  elevation  at  fixed  axial  coordinates  with  a  steepness  of  0.33  are 
presented  in  figure  8.  In  figure  9,  analytical  and  computed  results  for  a  wave  with  a  depth  of  2.01 
m  and  with  a  steepness  of  0.08  are  also  in  good  agreement. 


Figure  7.  Free-Surface  Elevations  at  Time  =  3.3  s  with  the  Sommerfeld 
Radiation  Condition  Used  at  the  Downstream  Boundary 


28 


Figure  8.  Free-Surface  Elevations  for  a  Steepness  of  0.33  at  (a)  x  =  0.44  m. 
(b)  x  =  4.87  m,  and  (c)  x  =  9.4  m  from  the  Upstream  Boundary  with  the 
Sommerfeld  Radiation  Condition  Used  at  the  Downstream  Boundary 


Figure  9.  Free-Surface  Elevations  for  a  Steepness  of  0.08  at  (a)  x  =  0.44  m, 
(b)  x  -  4.87  m,  and  (c)  x  =  9.4  m  from  the  Upstream  Boundary  with  the 
Sommerfeld  Radiation  Condition  Used  at  the  Downstream  Boundary 
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In  figure  10,  the  wave  was  started  with  the  fluid  at  rest  (</>  =  77  =  0  on  y  =  77 (x,r)).  A 
wavemaker,  programmed  as  in  equation  (50)  to  produce  an  Airy  wave,  is  placed  at  the  upstream 
boundary,  and  the  wave  is  allowed  to  propagate  freely  downstream.  The  free  surface,  initially  at 
rest,  is  only  developed  by  the  wavemaker.  The  downstream  boundary  is  four  wavelengths  from 
the  upstream  boundary.  The  water  depth  is  4  m  and  the  applied  upstream  Airy  wave  velocity 
potential  has  a  steepness  of  0.33.  These  are  the  same  parameters  as  were  used  in  figures  5  and  7. 
The  results  show  that  once  the  wave  is  developed  from  the  fluid  at  rest,  the  boundary  integral 
method  is  able  to  simulate  the  wave  over  a  long  period  of  time. 

In  figure  10,  the  wavemaker  was  positioned  upstream  throughout  the  computational  period. 
To  avoid  a  velocity  potential  jump  at  the  upstream/free-surface  interface  at  the  beginning  of  the 
experiment,  the  wavemaker  was  started  from  zero.  As  can  be  seen,  it  took  the  wavemaker  one 
wave  period  to  reach  the  steady  operating  condition  as  defined  by  the  Airy  wave.  With  time,  the 
wave  excitation  on  the  left  boundary  produces  a  wave  with  a  steepness  of  about  0.33  and  a  period 
of  1.26  s.  These  are  the  same  ffee-surface  wave  properties  obtained  from  linear  theory. 
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Figure  10.  Free-Surfa.ee  Elevations  for  a  Wave  Started  from  Rest  at  (a)  x  = 
0.44  m,  (b)  x  =  4.87  m,  and  (c)  x  =  9.4  m  from  the  Upstream  Boundary  with 
the  Sommerfeld  Radiation  Condition  Used  at  the  Downstream  Boundary 


CONCLUSIONS 


A  numerical  method  for  simulating  the  propagation  of  steep,  transient,  linear,  free-surface 
water  waves  has  been  presented.  The  method  is  based  on  boundary  integral  relations  utilizing 
linear  panels,  with  time  marching  performed  via  the  fourth-order  Adams-Bashforth  algorithm. 

This  algorithm  is  preferred  because  at  every  time  level  only  one  solution  to  the  Laplace  equation  is 
required.  To  obtain  starting  values,  a  predictor-corrector  method  was  used  during  the  first  three 
time  levels.  The  results  demonstrate  that  even  with  the  linear  panels  and  the  fluid  started  from  rest, 
the  numerical  method  is  capable  of  tracking  the  progression  of  steep,  free-surface  waves.  Also 
observed  is  a  close  agreement  of  these  results  with  the  analytical  linear  solution. 

It  is  important  to  note  that  the  proposed  method  does  not  depend  on  water  depth  or  on  wave 
properties  such  as  speed  or  wavelength.  This  capability  for  tracking  waves  generated  from  rest 
without  restrictions  on  wave  properties  will  be  essential  in  extending  the  method  to  the  fully 
nonlinear  water  wave  problem. 
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